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1  Introduction 


In  the  past  decades  several  different  long-range  imaging  systems  were  developed  working  in  both  the  visible 
and  infrared  spectral  bands.  A  direct  consequence  of  imaging  distant  scenes  is  the  effects  of  the  atmosphere. 
Especially  the  presence  of  turbulence,  which  become  non-negligible  and  affect  the  final  resolution,  limiting 
the  efficiency  of  identifying  a  possible  target.  Several  image  processing  algorithms  were  proposed  to  reduce 
the  impact  of  the  atmosphere  and  to  reconstruct  an  image  with  better  resolution  and  without  geometrical 
deformations,  see  [1,  15,  16,  19,  23,  21,  22,  24,  25,  28,  30,  32,  36,  42,  44,  45]  just  to  cite  a  few.  The  purpose 
of  this  project  was  to  explore  new  restoration  approaches;  to  improve  some  existing  techniques  (like  the 
Fried  deconvolution  [19]),  and  investigate  how  moving  target  detection  and  tracking  can  be  done  within  the 
framework  of  atmospheric  turbulence. 

2  Accomplishments 

During  this  two  years  project,  we  realized  three  main  tasks: 

•  We  created  a  open  dataset  of  static  and  dynamic  sequences  of  observations  through  atmospheric 
turbulence  in  order  for  the  community  to  create  fair  assessments  of  developed  algorithms.  An  un¬ 
dergraduate  student  was  involved  in  this  project.  This  work  was  published  in  J.  Gilles,  N.B.  Fer- 
rante,  "Open  Turbulent  Image  Set  (OTIS)",  Pattern  Recognition  Fetters,  Vol.86,  38-41,  2017,  DOI: 

10 . 1016/ j .patrec.2016. 12 . 020. 

•  Based  on  the  recent  idea  of  Fourier  Burst  Accumulation  deblurring,  we  proposed  several  combinations 
of  Wavelet  Burst  Accumulation  deblurring.  We  also  adapted  these  techniques  to  take  into  account  the 
non-rigid  deformations  due  to  the  turbulence.  This  work  was  published  in  J.  Gilles,  S.  Osher,  "Wavelet 
burst  accumulation  for  turbulence  mitigation".  Journal  of  Electronic  Imaging,  Vol.25,  No. 3,  033003- 
1-033003-9,  May  2016.  DOI:  10. 1117/1.  JEI.  25. 3. 033003. 

•  We  proposed  a  new  approach  to  detect  moving  objects  through  atmospheric  turbulence.  The  key  idea 
is  to  consider  geometric  and  oscillating  patterns  in  the  space+time  space.  This  work  involved  a  group 
of  undergraduate  and  graduate  students  during  the  Summer  2016  within  the  context  of  a  REU  (Re¬ 
search  Experience  for  Undergraduate)  project  funded  by  the  NSF  (my  colleague  Prof.  V.Ponomarenko 
is  the  program  director  in  our  department);  and  was  submitted  to  SIAM  Imaging  Sciences  Journal. 

The  following  sections  give  a  description  of  about  these  different  contributions. 

3  Open  Turbulence  Image  Set  (OTIS) 

Most  turbulence  mitigation  papers  in  the  literature  use  data  which  are  not  freely  available  hence  making 
fair  comparisons  difficult.  It  becomes  necessary  to  build  an  open  dataset  of  images  as  well  as  to  select 
some  metric  which  can  be  used  by  the  community  in  order  to  get  an  objective  comparison  of  the  different 
developed  algorithms.  In  this  section,  we  describe  such  an  open  dataset  we  created  in  San  Diego  thanks 
to  funds  from  a  university  grant  which  allowed  us  to  buy  the  basic  equipments.  This  dataset  is  made  of 
a  collection  of  sequences  of  static  scenes  as  well  as  dynamic  scenes  (i.e  with  a  moving  target).  Most  of 
the  static  sequences  correspond  to  the  observation  of  printed  charts.  Such  approach  permits  to  create  a 
groundtruth  associated  to  each  sequence  and  then  can  be  used  by  some  metric  to  assess  the  reconstruction 
efficiency.  We  categorized  the  observed  turbulence  as  either  “weak”,  “medium”  or  “strong”. 

The  work  described  in  this  section  was  published  in  J.  Gilles,  N.B.  Ferrante,  "Open  Turbulent  Image  Set 
(OTIS)",  Pattern  Recognition  Fetters,  Vol.86,  38-41,  2017,  DOI:  10.1016/j.  patrec  .2016.12.020. 
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3.1  Equipment  and  acquisition  procedure 

3.1.1  Equipment 

All  sequences  were  acquired  with  a  GoPro  Hero  4  Black  camera  modified  with  a  RibCage  Air  chassis  per¬ 
mitting  to  adapt  several  type  of  lenses.  We  always  used  a  25mm,  f/2.0  14d  HFOV  3MP  lens.  The  camera 
was  setup  at  a  1080p  resolution  and  a  framerate  of  24  frames  per  second  (fps).  The  camera  was  controlled 
by  the  standard  GoPro  App  on  a  Samsung  Galaxy  tablet. 

The  acquired  sequences  contain  both  natural  elements  from  the  observed  scene  as  well  as  artificial  "targets”. 
For  the  static  sequences,  we  used  two  charts  containing  some  geometric  patterns  at  different  spatial  frequen¬ 
cies  and  orientations  (see  left  two  images  of  Figure  1).  These  charts  were  printed  on  a  poster  (each  chart 
has  a  size  of  35x35cm)  and  held  up  by  a  homemade  wooden  stand.  For  the  dynamic  sequences,  we  used  a 
standard  remote  controlled  car  (see  right  image  of  Figure  1). 


Figure  1:  Two  left  images:  the  two  charts  serving  as  our  static  artificial  targets  after  being  printed  on  a 
poster.  Right  image:  remote  Controlled  car  utilized  as  our  moving  target  for  the  dynamic  sequences 


3.1.2  Procedures 

All  acquisitions  were  made  on  hot  sunny  days  in  order  to  guaranty  a  certain  amount  of  atmospheric  tur¬ 
bulence.  All  equipments  were  setup  on  a  practice  field  at  the  San  Diego  State  University.  This  field  is 
equipped  with  artificial  turf  which  reflects  very  well  the  sun  heat,  leading  to  high  level  of  turbulence.  The 
camera  stood  at  about  10cm  above  the  ground  observing  the  target  positioned  at  several  distances. 

After  all  acquisitions  were  done,  the  different  recorded  movies  were  downloaded  on  a  Linux  computer  and 
split  into  sequences  of  PNG  image  files  using  th ejfmpeg  command1.  The  different  region  of  interest  were 
finally  cropped  via  the  convert  command  (from  the  imagemagick  library2)  and  saved  as  individual  PNG 
sequences.  Since  the  Matlab®  software  is  widely  used  by  the  community,  we  also  provide  each  sequence 
as  a  Matlab  3D  matrix,  saved  in  a  .mat  file. 

Since  the  purpose  of  this  dataset  is  to  be  used  for  evaluating  turbulence  mitigation  algorithms,  all  sequences 
containing  the  two  above  mentioned  charts  arc  provided  with  a  groundtruth  image.  This  groundtruth  image 
contains  the  pristine  chart  after  being  downsampled  and  registered  to  the  actual  sequence.  In  practice,  we 
manually  registered  the  pristine  chart  on  a  temporal  average  of  the  input  sequence  using  the  GIMP3  software 
(this  procedure  is  summarized  in  Figure  2).  The  dynamic  sequences  arc  also  provided  with  their  respective 

'https : // f fmpeg . org/ 

2http : / /www . imagemagick . org/ 

3https : / / www . gimp .org/ 
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Figure  2:  Groundtruth  images  creation  procedure. 


groundtruths.  Each  groundtruth  is  a  sequence  of  binary  images  containing  the  bounding  box  corresponding 
to  the  moving  target  position.  These  groundtruths  were  created  with  the  software  Sensarea 4. 

3.2  Collected  data 

The  different  available  sequences  were  acquired  between  June  18th  and  August  16th,  2016.  We  pro¬ 
pose  OTIS  both  in  color  and  grayscales,  it  can  be  downloaded  for  free  at  https://zenodo.org/ 
communities/ otis/.  Characteristics  of  static  and  dynamic  sequences  are  given  in  Table  2  and  1, 
respectively.  Samples  frames  are  given  in  Figures  3  and  4. 


Figure  3:  Samples  of  frames  from  different  fixed  pattern  sequences  and  their  corresponding  groundtruths. 
The  first  row  corresponds  to  a  weak  turbulence  case  while  the  second  one  corresponds  to  a  strong  turbulence 
case. 

4http : / /www . gipsa-lab . grenoble-inp .fr/ -pascal. bertolino /bin/ download . php?f ile= 
sensarea . exe 
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Table  1:  Summary  of  the  different  dynamic  sequences 


Folder  Name 

Sequence 

Name 

Number 

of  im¬ 
ages 

Image  size 

Turbulence 

level 

Ground 

Truth 

Moving  Target 

Carl 

100 

200x200 

Medium 

Yes 

Car2 

315 

500x200 

Medium 

Yes 

Car  3 

51 

300x300 

Medium 

Yes 

Car4 

101 

300x300 

Medium 

Yes 

Table  2:  Summary  of  the  different  static  sequences 


Folder  Name 

Sequence 

Name 

Number 
of  im¬ 
ages 

Image  size 

Turbulence 

level 

Ground 

Truth 

Fixed  Background 

Door 

300 

520x520 

Strong 

No 

Fixed  Patterns 

Pattern  1 

64 

302x309 

Weak 

Yes 

Pattern2 

64 

291x287 

Weak 

Yes 

Pattern3 

300 

113x117 

Strong 

Yes 

Pattern4 

300 

109x113 

Strong 

Yes 

Pattern5 

300 

113x117 

Strong 

Yes 

Pattern6 

300 

109x113 

Strong 

Yes 

Pattern7 

300 

122x125 

Strong 

Yes 

Pattern8 

300 

119x122 

Strong 

Yes 

Pattern9 

300 

152x157 

Medium 

Yes 

Pattern  10 

300 

149x149 

Medium 

Yes 

Pattern  11 

300 

172x183 

Medium 

Yes 

Pattern  12 

300 

172x173 

Medium 

Yes 

Pattern  13 

300 

202x202 

Weak 

Yes 

Pattern  14 

300 

196x193 

Weak 

Yes 

Pattern  15 

300 

134x139 

Strong 

Yes 

Pattern  16 

300 

135x135 

Strong 

Yes 

Figure  4:  Frames  16  and  34  from  a  dynamic  sequence  and  their  corresponding  groundtruth  frames. 
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3.3  Quality  metric 


The  main  intent  of  OTIS  is  to  facilitate  the  comparison  of  turbulence  mitigation  algorithms.  Having  test 
sequences  as  well  as  their  groundtruth  images  is  not  sufficient  to  design  a  complete  evaluation  process. 
Indeed,  it  remains  to  choose  some  metrics  in  order  to  obtain  an  objective  comparison.  Denoting  In  and 
Iqt  the  restored  and  groundtruth  images,  respectively,  the  most  used  metric  to  compare  images  is  the  Mean 
Square  Error  (MSE)  defined  by 


MSE(Ir,  IGt)  =  ~^gt(®))2, 

X 

where  N  is  the  number  of  pixel  in  the  image.  Despite  its  popularity,  the  MSE  has  several  drawbacks,  notably 
it  doesn’t  take  into  account  the  geometric  information  contained  within  the  image.  The  Structural  Similarity 
Index  Measure  (SSIM)  was  proposed  by  [41]  to  circumvent  these  issues.  The  SSIM  is  a  perception-based 
model  which  consider  structural  (geometrical)  distortions  in  the  image.  Given  the  fact  that  one  of  the  major 
degradation  due  to  the  turbulence  is  the  geometrical  distortions,  the  SSIM  metric  appears  to  be  the  most 
adapted  metric  to  compare  turbulence  mitigation  algorithms. 

Regarding  the  assessment  of  tracking  algorithms  applied  on  the  dynamic  sequences,  the  methodology  de¬ 
veloped  in  [18]  can  be  used. 

4  Wavelet  burst  accumulation 

In  this  section,  we  describe  our  investigation  of  the  recently  proposed  weighted  Fourier  Burst  Accumulation 
(FBA)  method  into  the  wavelet  domain.  The  purpose  of  FBA  is  to  reconstruct  a  clean  and  sharp  image 
from  a  sequence  of  blurred  frames.  This  concept  lies  in  the  construction  of  weights  to  amplify  dominant 
frequencies  in  the  Fourier  spectrum  of  each  frame.  The  reconstructed  image  is  then  obtained  by  taking  the 
inverse  Fourier  transform  of  the  average  of  all  processed  spectra.  In  this  work,  we  first  suggested  to  replace 
the  rigid  registration  step  used  in  the  original  algorithm  by  a  non-rigid  registration  in  order  to  be  able  to 
process  sequences  acquired  through  atmospheric  turbulence.  Second,  we  proposed  to  work  in  a  wavelet 
domain  instead  of  the  Fourier  one.  This  led  us  to  the  construction  of  two  types  of  algorithms.  Finally,  we 
proposed  an  alternative  approach  to  replace  the  weighting  idea  by  an  approach  promoting  the  sparsity  in  the 
used  space. 

This  work  was  published  in:  J.  Gilles,  S.  Osher,  "Wavelet  burst  accumulation  for  turbulence  mitigation". 
Journal  of  Electronic  Imaging,  Vol.25,  No. 3,  033003-1-033003-9,  May  2016.  DOI:  10.1117/1.  JEI . 

25.3.033003. 

4.1  Fourier  burst  accumulation 

We  first  recall  the  concept  of  weighted  Fourier  Burst  Accumulation  (FBA)  proposed  in  [12]  and  we  introduce 
some  notations  which  will  be  used  throughout  this  section.  The  aim  of  the  FBA  method  is  to  retrieve  a 
deblurred  image  of  an  original  scene,  denoted  u,  from  a  sequence  of  observations,  assuming  that  those 
observations  arc  affected  by  hand-shake  blur.  We  will  denote  {vr  \^L ,  the  sequence  of  M  observations  such 
that,  Vi  =  1, . . . ,  M, 

Vi(x)  =  ( ki  *  m)(x)  +  n*(x),  (1) 

where  *  denotes  the  convolution  product,  x  the  position  in  the  image,  n;  some  noise  and  ki  is  a  blurring 
kernel  affecting  u  at  the  frame  i. 

The  authors  of  [12]  defined  the  weighted  Fourier  Burst  Accumulation  (FBA)  algorithm  in  the  following 
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way:  let  p  be  a  positive  integer,  the  restored  image,  up,  is  obtained  by 


'  M 


>(x)  =  FBAdviix)}^}^)  =  T  1  Wi{0F(vi)(0  (x), 


U=i 


Gff(|^(ni)(e)|p) 

with  wAt)  =  — - , 


(2) 


where  T  denotes  the  Fourier  transform  (£  are  the  frequencies)  and  Ga  is  a  Gaussian  filter  of  standard  devia¬ 
tion  a.  The  authors  mention  that  the  algorithm  is  not  sensitive  to  the  choice  of  a  and  it  can  be  automatically 
set  to  a  =  rn/J/oO  (where  mw  and  are  respectively  the  width  and  height  of  the  image).  In 

practice,  the  authors  apply  first  a  preprocessing  registration  step  (they  use  a  combination  of  SIFT  [26]  and 
ORSA  [14]  algorithms)  in  order  to  remove  affine  transformations  (translations,  rotations  and  homothety). 
The  authors  also  suggest  that  optional  denoising  and  sharpening  final  steps  can  be  applied  for  contrast  en¬ 
hancement  purposes. 

This  approach  is  based  on  the  following  principle:  Fourier  coefficients  which  represent  an  important  infor¬ 
mation  in  the  image  must  appeal-  consistently  in  all  observations.  Therefore,  the  weights  w,ip)  are  large  if 
Fourier  coefficients  F{vr)(p)  are  relatively  stable  through  time,  and  weak  otherwise.  Then  multiplying  the 
Fourier  spectrum  of  each  observation  by  these  weights  will  amplify  the  important  information  and  attenuate 
time  non-consistent  information.  Finally,  an  average  spectrum  is  computed  and  up  can  be  easily  retrieved 
by  inverse  Fourier  transform. 

As  mentioned  above,  the  authors  of  [12]  use  an  initial  rigid  registration  step  before  applying  the  FBA  al¬ 
gorithm.  Unfortunately,  for  atmospheric  turbulence  mitigation  purposes  the  observations  are  affected  by 
non-rigid  deformation  and  the  use  of  a  registration  step  based  on  SIFT  and  ORSA  algorithms  is  inefficient. 
In  [29],  a  non-rigid  regularization  technique  was  proposed  to  remove  atmospheric  non-rigid  distortions.  In 
particular,  the  deformation  fields  (I>,  can  be  estimated  via  some  optical-flow  algorithm.  We  kept  this  non- 
rigid  registration  step.  In  details,  it  consists  first  to  compute  an  average  frame  va(x)  =  2f=i  u(x)- 
Next,  the  deformation  mappings  <I>t,  such  that  Vi  =  1, . . . ,  M;  vt  =  4>,(na),  are  estimated  via  an  optical 
flow  algorithm.  Finally,  the  inverse  mappings  based  on  bilinear  interpolations  are  applied  to  each  frame: 
4>~ 1  fry)  — >  Vi.  This  registered  sequence  can  then  be  used  as  the  input  of  the  burst  accumulation  methods. 


4.2  Weighted  wavelet  burst  accumulation 

Other  image  representations  than  the  Fourier  representation  are  widely  used  in  image  processing,  notably 
wavelet  type  representations  (wavelets,  framelets,  curvelets,. . .).  Given  a  family  of  N  wavelets  {if)n}n=i> 
the  wavelet  representation  of  an  image  v  is  given  by  the  set  of  projections  {vn}^=1  =  {(v,  'ipn)\n-i  ■  In  this 
paper,  we  will  denote  W  a  wavelet  type  transform  of  the  image  v,  i.e  =  W(v).  The  inverse  wavelet 

transform  will  be  denoted  VV  ~ 1 ,  i.e  v  =  yV’_1({nn}^r=1). 

We  can  define  two  types  of  burst  accumulation  approaches.  The  first  one  processes  the  burst  accumulation 
directly  in  the  wavelet  domain,  we  will  denote  it  WWBA  (Weighted  Wavelet  Burst  Accumulation).  Here, 
the  weights  are  computed  from  the  wavelet  coefficients  themselves,  in  each  subband,  amplifying  the  most 
dominant  coefficients  through  time.  The  corresponding  formulation  is  given  by  equations  (3). 


M 


Up{x)  =  WWBA({Vi},  {V>n},p)  =  w  1  J  wj*(x)uj*(x) 


A=1 


N 


72=1  > 


(3) 


with 


ic?(x)  = 


GCT(K(x)n 
Ejii  GUKWb 
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The  Gaussian  filter,  Ga,  is  chosen  as  in  the  Fourier  case. 

The  second  type  of  accumulation  consists  in  doing  first  the  wavelet  decomposition  of  each  input  frame.  Then 
we  apply  the  Fourier  Burst  Accumulation  algorithm  from  [12]  (recalled  in  Section  4.1)  to  each  subband 
sequences  We  call  this  approach  the  Weighted  Wavelet  Fourier  Burst  Accumulation  (WWFBA) 

and  can  be  formulated  by  equations  (4). 

«p(x)  =  WWFBA({Vi},{i>n},p)  =  yW1  ({FBA({v?(x)}?ii,P)}n=i)  ■  (4) 

This  option  is  equivalent  to  deblur  the  wavelet  coefficients  themselves. 


4.3  Mathematical  characterization 

In  [12],  the  authors  characterized  the  FBA  algorithm  with  the  following  proposition. 


Proposition  1  Applying  the  FBA  algorithm  on  a  sequence 
to 

Up  =  ke*  u  +  n  where  ke  =  F~l 
and  n  is  the  weigh  ted  average  of  the  input  noise. 


provides  a  restored  image  up  equivalent 
/  M  \ 


^  WiF(ki )  , 


(5) 


u=i 


We  investigated  the  same  kind  of  characterization  for  the  two  previously  introduced  algorithms  in  the  wavelet 
domain.  First,  we  need  to  switch  to  the  filter  bank  point  of  view  about  the  wavelet  decomposition.  In  the 
previous  section,  we  denoted  the  wavelet  decomposition  of  v  by  {vn}^=1  =  {(v,  It  is  also  well- 

known  in  wavelet  theory,  using  the  frame  formalism  [10],  that  it  can  be  written  as  a  convolution  product: 


(6) 


where  ipn{x)  =  Un(— x).  The  inverse  wavelet  transform  is  given  by 

N 

v(x)  =  y\l1({vn}^=1)(x)  =  2  ( Vn  *  V’n)(x).  (7) 

71=1 


Therefore,  in  the  same  formalism,  the  wavelet  transform  of  each  observed  frame  vt  can  be  written  as  (we 
omit  x  in  the  following  in  order  to  simplify  the  notations) 


K"}n=l  =  iVi  *  VUn= 1  =  {(h  *U  +  m )  *  VUn=l  =  {h  *  U  *  fn}"=l  +  {tli  *  l/>n}%=1. 


\N  _ 


\N  .  = 


XN 


\  N 


(8) 


We  thus  obtain  the  following  characterizations. 


Analysis  of  the  WWBA  algorithm 

The  WWBA  restoration  algorithm  proposed  in  the  previous  section  is  equivalent  to 

M  N  M  N 

UP=SLlYj  (Wi(ki  *U*  ^n))  *  V’n  +  2  2  (^(n*  *  V’J)  *  Vw  (9) 

1=1  71=1  1=1  71=1 

Unfortunately,  pointwise  multiplication  and  convolution  do  not  commute  so  it  is  not  possible  to  write  this 
expression  as  the  convolution  of  u  with  some  equivalent  kernel  as  the  authors  did  in  [12]  for  the  Fourier 
burst  case. 
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Analysis  of  the  WWFBA  algorithm 

In  the  WWFBA  algorithm  case,  we  have  the  following  proposition. 

Proposition  2  Applying  independent  Fourier  burst  accumulations  on  each  subband  (n)  sequences  {v7jl}f£1 
is  equivalent  to 


N  M 

up  =  (ke  *  u)  +  n  where  =  Xj  X!  *  V’n  *  Vw  (10) 

n=  1  i=  1 


4.4  Non-linear  burst  accumulation 


The  frame  accumulation  techniques  described  in  section  4.2  are  all  based  on  linear  combination  (either 
by  pointwise  multiplications  or  convolutions)  of  each  frames  in  different  representation  domains.  These 
weights  basically  correspond  to  amplify  dominant  coefficients  thus  we  can  ask  if  it  is  possible  to  exploit  other 
approaches  to  perform  such  amplification?  We  proposed  an  alternative  to  the  use  of  weights.  Amplifying 
only  the  dominant  coefficients  can  be  interpreted  in  the  sense  that  only  those  dominant  coefficients  are 
important  to  represent  the  image.  Therefore,  instead  of  using  some  amplification,  we  can  imagine  to  keep 
the  dominant  coefficients  and  remove  the  other  ones.  In  other  words  the  restored  image  is  expected  to  have 
a  sparse  representation  in  the  used  representation  domains.  It  is  then  natural  to  promote  the  sparsity  in 
each  representation  before  doing  the  accumulation.  In  the  last  decade,  the  compressive  sensing  community 
widely  developed  such  concepts.  It  is  well  established  that  minimizing  models  based  on  L1-norm  provide 
sparse  solutions  [20].  Let  a  function  /,  the  simplest  model  to  find  a  sparse  representation  g  from  /  is  given 
by  (ll), 


g  =  argmin  ||^||i 
9 


2A II 9  * 


(11) 


It  is  well-known  that  the  solution  of  this  minimization  problem  is  given  by  soft-thresholding  /  with  threshold 
A,  and  is  given  by  (the  operators  are  understood  pointwise) 


g  =  Soft(f,  A)  =  X  max(|/|  -  A,  0). 


(12) 


Therefore,  denoting  V  the  representation  domain  which  can  be  either  T  or  VV,  we  propose  the  following 
general  sparse  burst  accumulation  model. 


rrp(x)  =  SVBAUviW}^},  A)  =  D”1  ^  Soft(V(Vi ),  A)  j  (x).  (13) 

Notice  that  in  the  case  of  a  wavelet-type  decomposition,  the  soft-thresholding  operator  is  applied  in  each 
subband. 

4.5  Experiments 

In  our  experiments,  we  decided  to  use  two  popular  families  of  wavelets  in  image  processing:  framelets 
and  curvelets.  Framelets  arc  basically  constructed  in  the  same  philosophy  as  for  classic  2D  tensor  wavelets 
except  that  no  downsampling  is  involved  in  the  process.  This  particularity  permits  to  guarantee  translation 
invariant  decompositions  which  is  important  for  image  processing.  Moreover,  it  is  well  known  that  such 
family  form  a  tight  frame  which  ensures  easy  reconstructions  and  permits  to  have  algorithms  which  arc  less 
sensitive  to  important  loss  of  information  (see  [35]  for  more  details).  Curvelets  also  form  a  tight  frame  but 
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their  main  particularity  is  in  the  fact  that  they  also  capture  directional  information  hence  providing  better 
representations  of  geometric  structures  in  images  (see  [5,  8]  for  details). 

In  order  to  distinguish  between  the  different  algorithms,  we  will  use  the  following  acronyms: 

•  FBA:  Fourier  Burst  Accumulation, 

•  Fr-WWBA:  Framelet  based  Weighted  Wavelet  Burst  Accumulation, 

•  C-WWBA:  Curvelet  based  Weighted  Wavelet  Burst  Accumulation, 

•  Fr-WWFBA:  Framelet  based  Weighted  Wavelet  Fourier  Burst  Accumulation, 

•  C-WWFBA:  Curvelet  based  Weighted  Wavelet  Fourier  Burst  Accumulation, 

•  SFBA:  Sparse  Fourier  Burst  Accumulation, 

•  Fr-SWBA:  Framelet  based  S parse  Wavelet  Burst  Accumulation, 

•  C-SWBA:  Curvelet  based  Sparse  Wavelet  Burst  Accumulation. 

We  experiment  the  proposed  algorithms  on  a  sequence  denoted  Barchart  1  (frames  1  and  25  arc  illustrated  in 
Figure  5).  We  ran  all  algorithms  (FBA,  Fr-WWBA,  C-WWBA,  Fr-WWFBA,  C-WWFBA,  SFBA,  Fr-SWBA 
and  C-SWBA)  without  any  registration  and  with  non-rigid  registration.  As  suggested  by  the  authors  of  [12], 
we  fix  p  =  11  for  all  weighted  methods.  For  the  sparse  based  methods,  the  parameter  A  is  set  to  A  =  0.5  for 
SFBA  and  A  =  0.001  for  Fr-SWBA  and  C-SWBA.  We  used  50  frames  in  each  test. 

The  Figure  6  presents  the  results  obtained  without  any  registration  step.  If,  compared  to  the  original  frames, 
cleaner  images  are  obtained,  it  is  worth  to  notice  that  some  geometric  distortions  remain  in  the  reconstructed 
images  thus  comforting  the  idea  that  a  non-rigid  registration  step  is  needed  to  deal  with  the  turbulence  prob¬ 
lem. 

The  results  obtained  when  using  the  non-rigid  registration  step  arc  given  in  Figure  7.  It  is  clear  that  the 
non-rigid  registration  step  is  essential  to  correct  the  geometric  distortions  induced  by  the  turbulence.  For 
instance,  we  can  notice  that  vertical  bars  in  the  barchart  arc  straighter  and  sharper  than  in  either  the  orig¬ 
inal  and  restored  without  non-rigid  registration  images.  Comparing  the  different  methods  with  non-rigid 
registration,  the  wavelet  based  options  give  clearer  restored  images  than  the  FBA  technique.  The  advantage 
looks  like  to  be  in  favor  of  the  sparse  options  as  they  provide  images  with  more  contrast  than  the  weighted 
approaches. 


*1  Mill  I  %  llii* 


Figure  5:  Original  sequence  Barchart  1.  The  consecutive  columns  correspond  to  the  frames  1  and  25, 
respectively. 
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FBA 


Fr-WWBA 


Figure  6:  Restoration  results  on  Barchart  1  sequence  (50  frames)  without  non-rigid  registration  and  p  =  11 
(without  any  final  denoising  or  sharpening  step). 
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FBA  Fr-WWBA 


C-WWBA  Fr-WWFBA 


C-WWFBA  SFBA 


Fr-SWBA  C-SWBA 


Figure  7:  Restoration  results  on  Barchart  1  sequence  (50  frames)  with  non-rigid  registration  and  p  =  11 
(without  any  final  denoising  or  sharpening  step). 
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5  Moving  target  detection  through  turbulent  media 


In  this  section,  we  investigate  how  moving  objects  can  be  detected  when  images  are  impacted  by  atmo¬ 
spheric  turbulence.  We  present  a  geometric  spatio-temporal  point  of  view  of  the  problem  and  show  that  it 
is  possible  to  distinguish  movement  due  to  the  turbulence  vs.  moving  objects.  To  perform  this  task,  we 
propose  to  extend  2D  cartoon+texture  decomposition  algorithms  to  3D  vector  fields. 

This  work  has  been  submitted  to  SIAM  Imaging  Science  Journal:  J.  Gilles,  F.  Alvarez,  N.  Ferrante,  M. 
Fortman,  L.  Tahir,  A.  Tarter,  A.  von  Seeger,  “Detection  of  moving  objects  through  turbulent  media.  Decom¬ 
position  of  Oscillatory  vs  Non-Oscillatory  spatio-temporal  vector  fields.”. 

5.1  Moving  objects  in  turbulent  media 


Figure  8:  Several  frames  from  the  OTIS  Carl  sequence.  The  second  row  shows  the  velocity  vector  field 
corresponding  to  each  frame. 

We  start  to  illustrate  the  specificity  of  sequences  impacted  by  the  presence  of  atmospheric  turbulence. 
Figures  8  and  9  present  several  frames  from  the  sequences  Carl  and  Car2  from  the  Open  Turbulent  Image 
Set  (OTIS)  [17]  described  in  Section  3.  Clearly,  the  turbulence  creates  dynamic  local  deformations  which 
are  confirmed  by  the  velocity  vector  fields  given  on  the  second  row  of  each  of  those  figures  (the  color  wheel 
on  the  left  gives  the  correspondence  between  the  color  and  the  vector  direction).  All  pixels  of  all  frames  are 
clearly  subject  to  some  movement.  Moreover,  these  vector  fields  are  mainly  made  of  random  vectors  (both  in 
direction  and  magnitude).  Therefore,  it  is  almost  impossible  to  distinguish  the  movement  information  which 
corresponds  to  moving  object  vs.  moving  atmosphere,  therefore  a  new  approach  is  needed.  The  standard 
vision  of  the  problem  is  biased  by  the  fact  that  we  watch  the  sequence  frame  by  frame  (i.e.  a  temporal 
succession  of  2D  images).  We  suggested  to  consider  the  whole  sequence  as  a  single  spatio-temporal  space. 
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Figure  9:  Several  frames  from  the  OTIS  Car2  sequence.  The  second  row  shows  the  velocity  vector  field 
corresponding  to  each  frame. 


In  particular,  if  we  visualize  the  colored  vector  field  as  a  3D  cube  of  data  (i.e.  each  pixel  in  the  3D  space  is 
a  vector),  we  can  see  the  emergence  of  geometric  patterns.  Indeed,  if  the  turbulence  movement  is  random 
in  both  space  and  time,  the  movement  of  a  moving  object  is  expected  to  be  “consistent”,  i.e  a  moving 
object  follows  some  trajectory  which  is  not  supposed  to  be  erratic.  This  concept  is  illustrated  in  Figure  10 
where  slices  corresponding  to  x  —  time  planes  (y  is  fixed)  are  given.  This  visualization  confirms  that 
a  moving  object  creates  a  geometric  pattern  in  the  space-time  volume  while  the  atmospheric  turbulence 
velocity  field  behaves  like  an  oscillating  function.  Therefore,  we  suggest  that  detecting  moving  objects 
through  atmospheric  turbulence  can  be  achieved  by  separating  non-oscillating  vs.  oscillating  components  of 
the  vector  fields.  This  idea  is  explored  in  the  next  section. 


(a)  Carl  (b)  Car2 


Figure  10:  3D  spatio-temporal  patterns  visualization  in  the  velocity  vector  fields  of  the  Carl  and  Car2 
sequences. 


5.2  Decomposition  methods 

The  idea  of  separating  oscillating  vs.  non-oscillating  components  has  been  widely  studied  in  image  process¬ 
ing.  The  specific  purpose  of  these  investigations  was  to  separate  the  cartoon  component  from  the  texture 
component  in  a  given  image.  The  cartoon  and  texture  parts  correspond  to  non-oscillating  and  oscillating 
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function,  respectively.  In  this  work,  we  extended  such  type  of  algorithm  to  spatio-temporal  vector  fields. 
We  first  recall  the  formalism  used  in  the  image  processing  case  and  next  we  introduce  the  vector  field  case. 


5.2.1  Image  decomposition 

Given  an  image  /,  the  purpose  of  the  cartoon  +  texture  decomposition  is  to  find  two  images  u  (the  cartoon, 
i.e.  piecewise  constant  function)  and  v  (the  textures,  i.e.  oscillating  function)  such  that  f  =  u  +  v.  This 
concept  was  proposed  by  Meyer  [31]  in  his  investigation  of  the  Rudin-Osher-Fatemi  (ROF)  model  [34]. 
Several  authors  [2,  3,  4,  33,  38,  40]  then  proposed  different  modifications  of  Meyer’s  model  leading  to 
tractable  algorithms.  In  particular,  the  case  of  using  Besov  spaces  [39],  Bsp  q  (we  consider  the  homogeneous 
version  of  these  spaces),  is  numerically  interesting  since  Besov  spaces  can  be  characterized  by  the  mean  of 
wavelet  coefficients  [9,  1 1,  27],  Using  the  fact  that  B\  x  c  BV  and  G  c  B^x  OD,  we  can  formulate  a  Besov 
based  variational  decomposition  model  by 


where 


(it,  v )  =  arg  min 


ueB[  1,veBa_ 


Bltl 


+  (2A)  1\\f  —  (u  +  v)\\L2, 


j* 


B 


if  Ill’ll  oco  n 

l,oo 

otherwise 


(14) 


(15) 


Such  model  can  easily  be  solved  by  iterating  soft  thresholding  of  the  wavelet  coefficients.  Denoting  i'V(f  ) 
the  wavelet  transform  of  /  (we  assume  a  real  transform)  and  Shrink(x,  A)  =  sign(x)  max(0,  x  —  A),  the 
numerical  algorithm  can  be  resumed  by  Algorithm  1 . 


Algorithm  1  Besov  based  decomposition  numerical  algorithm. 

Inputs:  image  to  decompose  /,  parameters  A,  //,  maximum  number  of  iterations  Nmax 
Initialization:  n  =  0,  rr°  =  n°  =  0 

repeat 

vn+1  =  f  _  un  _  yy-1  ( Shrink(W(f  —  un),  2/i)) 
un+ 1  =  W_1  (Shrink(W(f  -  nn+1),2A)) 
until  max  (|j'Un+1  —  un\\L2,  ||nn+1  —  || z,2 )  <  10~6  or  Nmax  iterations  are  reached 

return  un+1,vn+1 


5.2.2  Spatio-temporal  vector  field  decomposition 

Next,  we  extended  the  previous  decomposition  model  to  spatio-temporal  vector  fields.  Let  us  denote 
f(i,j,  n)  the  input  sequence  of  N  frames  where  i,j  arc  the  2D  spatial  variables  and  n  is  the  frame  number. 
The  decomposition  model  can  be  easily  extended  to  3D  by  simply  using  a  3D  wavelet  transform  instead  of 
the  2D  one,  while  the  algorithm  itself  remains  unchanged.  However,  we  want  to  decompose  vector  fields 
corresponding  to  the  estimated  velocity  from  frame  to  frame  hence  the  previous  algorithm  must  be  adapted 
to  vector  fields.  If  we  denote  the  spatio-temporal  vector  field  v(i,j,n)  =  (vi(i,j,n),V2(i,j,n)),  the  first 
idea  which  comes  to  mind  is  to  process  the  vector  field  components  v\  and  vo  separately  and  then  recompose 
the  final  vector  field.  Unfortunately,  this  approach  will  not  perform  well  in  some  cases,  like  when  one  of  the 
components  is  oscillating  while  the  other  is  constant.  We  need  to  process  the  vector  field  as  a  single  object. 
We  can  simply  achieve  this  by  rewriting  the  vector  field  as  a  complex  vector  field,  i.e 

v(i,j,n)  =  vi  (i,j,n)  +  iv2(i,j,n). 
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Now,  we  can  use  the  same  3D  decomposition  algorithm  only  if  we  change  the  standard  wavelet  transform 
by  a  complex  wavelet  transform  and  then  use  the  complex  version  of  the  thresholding  operator  defined  by 

\/z  =  \z\el6  e  C  ,  CShrink(z ,  A)  =  max(0,  \z\  —  X)eld . 


Moreover,  we  suggest  a  last  modification  to  the  algorithm  to  better  take  into  account  geometric  patterns 
with  specific  orientations.  Indeed,  as  shown  in  the  Figure  10,  the  pattern  associated  to  the  moving  object 
correspond  to  a  piecewise  constant  function  with  a  direction  not  parallel  to  the  space  axis.  Therefore,  it  is  of 
interest  to  use  a  representation  which  takes  into  account  such  different  orientations.  The  3D  classic  wavelet 
transform  is  built  on  a  tensor  approach,  i.e.  is  separable  into  three  ID  transforms  with  respect  to  each  axis, 
hence  does  not  allow  to  take  care  of  different  orientations.  An  alternative  to  the  classic  wavelet  transform  is 
the  curvelet  transform  [5,  6,  7,  8,  13,  37,  43]  where  the  curvelet  filters  arc  built  on  angular  wedges  in  the  3D 
Fourier  space.  Since,  curvelets  arc  wavelet  type  functions,  it  is  possible  to  formalize  curvelet  based  Besov 
spaces  which  we  will  denote  C®  .  We  can  then  replace  the  use  of  Besov  spaces  B^  q  by  curvelet  spaces  Cp.q 
and  the  curvelet  based  decomposition  model  consists  in  solving  model  (16). 


where 


(u,v) 


argmin  ||rt|jpi 
ueCl^veC^ 


+  (2A)  1||/  —  (u  +  v)\\L2, 


T* 

JC 


if  IMI^  <  H 

otherwise 


(16) 


(17) 


If  we  denote  C  the  curvelet  transform  of  /,  the  numerical  algorithm  which  solves  model  (16)  is  given  in 
Algorithm  2. 


Algorithm  2  Curvelet  based  decomposition  numerical  algorithm. 

Inputs:  image  to  decompose  /,  parameters  A,  //,  maximum  number  of  iterations  Nmax 
Initialization:  n  =  0,  rr°  =  v°  =  0 

repeat 

vn+l  _  y  _  un  _  q-1  (cShrink(C(f  —  u11),  2/x)) 
un+ 1  =  C~l  ( Shrink{C{f  -  vn+1),2X)) 

until  max  (|jrrn+1  —  un\\L2,  ||nn+1  —  vn\\L2Sj  <  10~6  or  Nmax  iterations  are  reached 
return  un+1,vn+1 


5.3  Experimental  results 

In  this  section,  we  present  the  results  obtained  by  Algorithm  2  on  the  sequences  Carl  and  Car2  from  the 
OTIS  dataset  presented  in  Section  3.  We  experimentally  observed  that  the  choice  A  =  n  =  1  works  well 
for  all  sequences  and  hence  was  used  in  all  our  experiments.  The  maximum  number  of  iterations  was  fixed 
to  Nmax  =  5.  Figures  11  and  12  provide  the  geometric  and  oscillating  components  of  the  velocity  vector 
fields  corresponding  to  each  sequence,  respectively.  As  expected,  the  geometric  component  highlights  the 
movement  of  the  moving  object  while  the  oscillating  component  contains  the  turbulence  movement.  Even 
on  the  difficult  sequence  Car2  where  the  object  is  small,  moves  close  to  an  edge  from  the  background 
and  hence  is  strongly  affected  by  the  turbulence,  the  algorithm  does  a  good  job.  Furthermore,  in  order  to 
have  a  comparison  with  Figure  10,  we  provide  the  spatio-temporal  plots  for  the  extracted  components  in 
Figure  13.  Again,  it  is  obvious  to  see  that  the  expected  geometric  patterns  corresponding  to  moving  objects 
are  enhanced  in  the  geometric  component  making  them  easier  to  detect. 
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Figure  11:  Decomposition  results  on  the  OTIS  Carl  sequence.  The  first  row  corresponds  to  the  geometric 
component  and  the  second  row  the  oscillating  component. 


9 

Figure  12:  Decomposition  results  on  the  OTIS  Car2  sequence.  The  first  row  corresponds  to  the  geometric 
component  and  the  second  row  the  oscillating  component. 


(a)  frame  50 


(b)  frame  150 


(c)  frame  200 
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(a)  Carl  -  Geometric  component 


(b)  Carl  -  Oscillating  component 


Figure  13:  3D  spatio-temporal  patterns  visualization  of  the  geometric  and  oscillating  components  of  the 
Carl  and  Car2  sequences  velocity  vector  fields. 


6  Work  in  progress 

We  arc  currently  making  progress  on  two  topics:  parameters  estimation  of  the  Fried  deconvolution  and  the 
design  of  a  lucky  imaging  super-resolution  algorithm. 

•  The  Fried  kernel  models  the  impact  (in  terms  of  blur)  of  the  atmospheric  turbulence  in  images  and 
was  used  in  [19]  to  perform  the  final  deconvolution  needed  to  get  sharper  images.  The  main  issue 
with  the  current  Fried  formulation  is  the  need  to  know  four  parameters.  We  started  some  investigation 
on  mathematically  simplifying  the  kernel  with  the  goal  of  building  an  efficient  method  to  estimate  the 
remaining  parameters.  Unfortunately,  our  first  experiments  showed  that  the  problem  is  much  more 
complex  than  expected.  In  particular,  priors  used  to  characterize  the  parameters  are  non-convex.  We 
recently  started  a  collaboration  with  Dr.  Charles  Deledalle  from  France  who  just  started  a  sabbatical 
year  at  UCSD.  Dr.  Deledalle  is  a  specialist  of  parameter  estimation.  We  are  starting  to  investigate 
new  approaches  to  solve  this  problem. 


18 

DISTRIBUTION  A:  Distribution  approved  for  public  release. 


•  The  lucky  imaging  idea  was  proposed  many  years  ago  but  in  our  opinion,  it  is  completely  under¬ 
exploited.  A  SDSU  graduate  student,  Francis  Alvarez,  is  currently  working  on  this  topic  for  his 
master’s  thesis.  Our  idea  is  that  the  atmosphere  locally  behaves  like  an  optical  lens.  The  presence  of 
turbulence  locally  modify  this  behavior  in  the  sense  that  the  atmosphere  can  correspond  to  either  a 
zoom  in  or  zoom  out  effect.  A  good  idea  is  then  to  keep  only  the  information  corresponding  to  a  zoom 
in,  while  getting  rid  of  the  information  from  zoom  out  regions.  Doing  so,  we  can  expect  to  keep  only 
the  best  information  from  the  sequence  to  reconstruct  a  restored  image  and  we  even  expect  to  be  able 
to  do  some  very  efficient  super-resolution.  As  of  today,  we  already  found  how  to  distinguish  zoom 
in  vs.  zoom  out  regions  in  a  given  frame.  We  arc  currently  working  on  how  this  information  can  be 
incorporated  in  a  super-resolution  algorithm. 
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